In this practical, we will analyse data pertaining to customer waiting times in a governmental office. The data provided to us is the \(waiting.csv\) file are daily average waiting times from January 3, 2017 to October 26, 2018.
Our goal is to model and monitor the time customers spend waiting, and to take action if the service quality deteriorates too badly, i.e. if waiting times are unacceptably high.
Read in the data and plot the daily waiting times. Compute and display the five-number summary (boxplot) of the daily waiting times for each weekday. What do you observe?
| Weekday | Min | Q1 | Mean | Q3 | Max |
|---|---|---|---|---|---|
| Friday | 148 | 342.25 | 485 | 622.00 | 918 |
| Monday | 116 | 297.00 | 446 | 558.75 | 999 |
| Thursday | 81 | 217.00 | 373 | 487.25 | 938 |
| Tuesday | 90 | 262.00 | 431 | 550.00 | 1154 |
| Wednesday | 95 | 293.25 | 471 | 597.00 | 1182 |
\(Conclusion\) : our main observations regarding the charts below are the following:
Using a normal approximation produce an upper control limit or confidence line for the waiting times at the one-year return level. What do you notice? Using a Q-Q plot, comment on whether the Normal model is appropriate.
########Assuming that waiting times are normally distributed
#95% CI--> alpha = 0.05 We can get z(alpha/2) = z(0.025) from R:
nrow<- nrow(data)
qnorm<- qnorm(.975) #1.959964
meanData<- mean(data$Average.Wait.Seconds) #440.6398
sdData<- sd(data$Average.Wait.Seconds) #206.19
###########Our margin of error is
upper_control_limit<- mean(data$Average.Wait.Seconds)+
qnorm*(sdData/sqrt(nrow))
#If X~N(μ, σ)--> pnorm(x, μ, σ) function to calculate P(X > upperbound).
# r=1
probability_exceed_UCL<- 1 - pnorm( upper_control_limit,
meanData, sdData)
probability_exceed_UCL
## [1] 0.4630697
\(Conclusion\): our main observations regarding the analisis are:
The aim of this analysis is to focus on negative effects (long wait times). Explain how you would proceed using 1) a block maxima and 2) a peaks-over-threshold approach. For each method, carry out the required data aggregation and transformation.
We have proposed the following steps:
Block Maxima: \[M_{i}=max\{X_{(i−1)n+1},\cdots,X_{in}\},i=1,⋯,m\]
Then, we fit the a Generalized extrem Value, assuming that it follows \(GEV(\mu_{n},\sigma_{n},\xi)\) distribution and we will focus on maximize the log-likelihood and the estimation of the parameters.
In this step, we decide to test if the parameters correspond to the Gumbel distribution.
Finally, we asses the model and compute m-year return levels.
#we should first extract the maximum temperature values--> block maxima.
#week variable creation
# derive AMS for maximum precipitation
data_ts <- xts::xts(data$Average.Wait.Seconds, order.by = data$Date) #time series creation
wms_data <- xts::apply.weekly(data_ts, max) #create the weekly maximal serie
#According to the Fisher–Tippett–Gnedenko theorem, the distribution of block maxima can be approximated by a generalized extreme value distribution.
gev.fit2 <- extRemes::fevd(as.vector(wms_data),method = "MLE", type = "GEV")
| Location | Scale | Shape | |
|---|---|---|---|
| Estimates | 585.94 | 170.64 | -0.13 |
| Std.errors | 19.87 | 14.06 | 0.07 |
######Testing the Gumbel distribution (ξ=0) hypothesis with LRT#####################
#################################################################################
gev_gumb <- fevd(as.vector(wms_data), type="Gumbel", units="deg C")
test_gumbel<- lr.test(gev_gumb,gev.fit2)
test_gumbel
##
## Likelihood-ratio Test
##
## data: as.vector(wms_data)as.vector(wms_data)
## Likelihood-ratio = 2.478, chi-square critical value = 3.8415, alpha =
## 0.0500, Degrees of Freedom = 1.0000, p-value = 0.1154
## alternative hypothesis: greater
#leads us to not reject the Gumbel hypothesis !
########################################################################
\(Conclusion\) : our main observations are the following:
we are going to perform the analysis of the plot of the fit in the question D.
We have proposed the following steps:
First, we need stablish a large enough treshold \(u\), where the following formule need to be fullfill: \[\mathbb{P}\left(X_{i}-u>y|X_{i}>u\right)\] In order to do that, we are going to run a mean residual plot and then we will select the value where the plot start looks linear.
Then, we fit the all the dataset to a Generaliaze Pareto distribution.
Finally, we asses the model and compute m-year return levels.
# maximum likelihood estimation
pot.fit <- fevd(as.vector(data_ts), method = "MLE",
type="GP",
threshold = threshold)
\(Conclusion\) : our main observations are the following.
we are going to perform the analysis of the plot of the fit in the question D.
Propose a model for the data you have processed in the previous question. Make sure to justify your model choice using residuals and goodness-of-fit tests. (Hint: you may use the evd package.)
Now, we will use likelihood-ratio test for two nested extreme value distribution models,
the function used comes from the extRemes package.
par(mfrow=c(1,2))
plot(gev.fit2, type = "qq2")
plot(pot.fit, type = "qq2")
Figure 1: Quantile plot GEV-GP
\(Conclusion\): our main observations regarding the analisis are the following.
Finally, use your model to derive an upper control limit or confidence line for the waiting times at the one-year return level.
| name | values |
|---|---|
| one year return level | 1140.82 |
| Lower bound | 1041.54 |
| Upper bound | 1240.09 |
\(Conclusion\) :
In this second practical, we will look at product sales data belonging to a large American retailer.
As the person in charge of re-stocking products after they are sold, you would like to avoid a ‘stock-out’, that is, when on a given day there are not enough products to satisfy customer demand.
Consider the observations recorded in sales 1.txt, which contain the daily sales volumes of Product 1 over n = 1913 days.
Construct a figure which shows the daily sales volume, and study the values just before the product is out of stock. What do you notice? Explain why this is related to stock-out.
plot(sales$Sales, type = "l", main = "Daily Sales")
outofstock <- which(diff(sales$Sales) == -sales$Sales & sales$Sales !=0)
## Warning in diff(sales$Sales) == -sales$Sales: longer object length is not a
## multiple of shorter object length
abline(v = outofstock, col = "red")
sales %>% dplyr::mutate(date = as.Date(1:1913)) %>% dplyr::filter(Sales == 0)
## Sales date
## 1 0 1970-11-28
## 2 0 1971-01-20
## 3 0 1971-01-21
## 4 0 1971-11-29
## 5 0 1973-06-29
## 6 0 1973-06-30
## 7 0 1973-11-28
## 8 0 1974-11-28
It seems there is some seasonality element in the data, as the product is out of stock alomsto every year at the same time (around the end of November for almost every year), hence it could be avoided as it is repeated thorugh different years.
Assume Product 1 is a perishable good sold at a constant price of 10.– over the time period considered. Items of Product 1 are produced at cost 1.– and items which are unsold can be salvaged at 10% of cost, that is, at value s = 10 ct. Suppose that the sales of Product 1 are distributed according to F. Using the newsvendor model, give a formula for the critical fractile q = F−1(p).
Otherwise, we can calculate it by hand.
#price
pr <- 10
#cost of production
c <- 1
#salvage cost
s <- 0.1
#cost of shortage
cs <- pr-c
#cost of overage
co <- c-s
#critical fractile
p <- cs / (cs + co)
In this case, the formula to calculate the critical fractile is the following:
\[ q = F^-1 (p)\]
is equal to
\[ q = F^{-1} (\frac {cs}{cs + co}) \]
The cs as the cost of shortage, which corresponds to the cost per unit it would incur if there is still demand but we are out of stock, hence a missed sale which correspond to the price minus the costs (p - c). On the other hand co is the cost of overage, which is the cost per unit of having too many products compared to the demand, hence a product produced and unsold which is the cost of production minus the salvage cost (c - s).
Explain how the critical fractile is related to the Value at Risk, and give an interpretation of the Expected Shortfall at the level p.
The critical fractile, balances the cost of being understocked and the total costs of being either overstocked or understocked, and it gives the optimal quantity that should be produced to minimizes the costs. It’s the quantile which optimize the probability that the demand is lower than the quantity Q that has been ordered.
The Value at Risk (VaR) denotes a quantile of the distribution, which is the smallest number such that the random variable X that we are studying exceeds a value x with probability lower than 1 - alpha. In other words, with probability alpha, the variable X will be smaller than VaR. Hence, we are looking at the heavy tails of the distribution of the variable.
The interest is when the newsvendor model gives a high quantile to look at (e.g. above 0.9, like in our case in the previous question). What we want to avoid is to go beyond the critical fractile, because in that case we will not choose the optimal quantity, and we will increase the costs and hence decrease the profit, even though this will happen with a low probability (of 1 - q). This is linkable to the Value at Risk as it is also describing an extreme situation in which the variable takes a value that is higher than a certain threshold, and, as we are in the risk analysis field we are studying the losses on investments, we want to avoid these extreme cases, as they mean a high loss.
The Expected Shortfall at level alpha is the average value of X when it is higher than VaR. It is related to VaR by
\[ ES_\alpha (x) = \frac{1}{1 - \alpha} \int_\alpha^{1}q_u (F_X)du = \frac{1}{1 - \alpha} \int_\alpha^{1} VaR_\beta(L)d\beta \]
The Expected Shortfall at level p will then be the value of the sales if they are higher then the VaR, so for the level p, the ES will be the average value of the sales that are higher than the level p that is calculated as the quantile of the critical fractile q.
Fit a Poisson model to the sales volume data by assuming that Yi ∼ Poisson(μ) and using maximum likelihood estimation for μ. (Hint: you may use MASS::fitdistr). According to this model, what is the estimated critical fractile?
poisson <- fitdistr(sales$Sales, "Poisson")
#lambda of the poisson distribution
lambda <- poisson$estimate
#quantile
q <- qpois(p, lambda)
In the case of a poisson distribution, we find a critical fractile by using the following formula:
\[ q = F^{-1} (p) = \sum_{i = 0}^p \frac {e{-\lambda}\lambda^i}{i!}\]
And it corresponds to 114.
Suppose we model the sales as an extreme value distribution using a peaks-over-thresholds method. Give the Value at Risk at level p under this model.
#Mean Residual Life Plot
extRemes::mrlplot(sales$Sales, main="Mean Residual Life Plot")
#Finding and appropriate threshold
extRemes::threshrange.plot(sales$Sales, r = c(100, 200), nint = 120)
#let's fix the threshold at 140 as after it the distribution seems to be flat in the second graph and there is the elbow in the first one
u <- 190
plot(sales$Sales, type = "l", main = "Daily Sales") +
abline(h = u, col = "red")
## integer(0)
#Weights for the observations above the threshold chosen
W <- sales %>% dplyr::filter(Sales > u) %>% dplyr::mutate(W = Sales - u) %>% dplyr::select(W)
#fitting a GPD, method amle as we want the shape parameter positive
gpd <- gPdtest::gpd.fit(W$W, method = "amle")
#extracting the parameters
#shape
beta <- gpd[1,]
#scale
xi <- gpd[2,]
#setting the alpha for the VaR to the level we have calculated
alpha <- p
#calculating Fbar with the approximation
Fbar <- nrow(W) / nrow(sales)
#calculating the values at risk
VaR <- (u + (beta/xi) * ( ( ( (1 - alpha) / (Fbar) ) ^ (-xi) ) - 1) )
We wanted to use the first two graphs to determine the therhsold that we want to use. What we can see in the first one is that between the value 100 and 150 the curve seems to flatten, which is confirmed by the second graphs, as the parameters are becoming constant already from before 100. However, choosing a threshold between these two values would give too many observations above it, in our opinion, hence we decide to move it a little bit further up. We decided to be safe and to choose a threshold of 190, so that we have a quantile that is around the 80%.
We move on by calculating the weights of the observations that are above this threshold (that are the points in the graph above the red line).
Then, using these observations we fit a GPD and we estimate the parameters, that we will insert in the formula to calculate the VaR.
The resulting value is 189.9589851.
Perform a binomial back-test over the last 300 days in the dataset for the models in d) and e), using a window size of 365 days. What are your conclusions?
To backtest: 1. Split the data and create a test set including the last 300 observations of the dataset 2. Calculate the VaR for each one of the observations inside the test set using a window of 365 days 3. Compare the estimated VaR to the real values and see if there is any violation, meaning that the real value is higher than the VaR for a given day 4. Consider that the violations follow a binomial distribution(300, 1 - p), test the hypothesis that the observed proportion of violations over the test set is different from the expected proportion given by the binomial distribution with a binomial test
We will split the data in two subgroups, the first one being the training set which we will use to estimate the VaR will comprehend the observations up to the last 300 days of the data, while the last 300 observations will be the test set to which we will compare the VaR we will find.
salesTest <- sales$Sales[(nrow(sales)-299):(nrow(sales))]
Then, we calculate the Value at Risk by fitting a poisson distribution on the window preceeding the observation in which we are interested.
We start by creating a function that works in the following way: 1. use the same probability as before (hence 0.91) 2. fit the poisson on the window 3. determines the lambda of the poisson distribution 4. caclulate the VaR by taking the 0.91 quantile of the poisson we fit
poissVaR <- function(df){
#definition of the probability
p <- 0.91
#fitting the poisson
poisson <- fitdistr(df$Sales, "Poisson")
#lambda of the poisson distribution
lambda <- poisson$estimate
#quantile, which corresponds to the VaR
VaRpoiss <- qpois(p, lambda)
}
Then, we can create a loop to get the VaR for the observations inside the test set.
#create a vector to store the values
poissVAR <- as.numeric(1:300)
for(i in 1:300){
#select window
wind <- sales %>% dplyr::slice((1913-i):(1913-i-365))
#calculate VaR
poissVAR[i] <- poissVaR(wind)
}
poissVAR
## [1] 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106
## [19] 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106
## [37] 106 106 106 106 106 106 106 106 106 106 106 106 106 106 105 105 105 105
## [55] 105 106 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105
## [73] 105 105 105 105 105 105 105 105 105 105 105 106 105 105 105 105 105 106
## [91] 106 106 106 105 105 105 106 106 106 105 105 105 105 106 106 106 105 105
## [109] 105 105 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106
## [127] 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 106 105 105
## [145] 106 106 106 106 106 106 105 105 105 106 105 105 105 105 105 105 106 106
## [163] 106 105 105 105 105 106 106 105 105 105 105 105 105 105 105 105 105 105
## [181] 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105
## [199] 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105 105
## [217] 105 105 105 105 105 105 105 105 105 105 105 106 106 106 106 106 106 106
## [235] 106 106 106 107 107 106 107 107 107 107 107 107 106 106 106 106 106 106
## [253] 106 106 106 106 106 106 107 107 107 107 107 107 107 107 107 107 107 107
## [271] 107 107 107 107 107 107 107 107 107 107 107 107 107 107 107 107 107 107
## [289] 107 107 107 107 107 107 107 107 107 107 107 107
As we can see they take they almost take the same values, this is due to the fact that the window is quite large (365), hence the change from one window to the subsequent one will not be very big, as we just drop the furthest away observation to gain the most recent one that is the day before the one we are considering. Hence, the fitting of the poisson distribution will be almost the same, as it will be the lambda and the quantile.
Now, we compare the VaR we estimated to the real values, to see how many violations there are. The should be following a bernoulli distribution with n = 300 and p = 1 - alpha.
#number of observed violations
poiss.violations <- sum(ifelse(poissVAR > salesTest, 0, 1))
#number of theoretical violations, given by the mean of the bernoully distribution
ev <- 300*(1 - p)
We find out that there are 89 violations.
We can compare this value to the theoretical one that we expect to get, which is 27.2727273.
We can test whether they are statistically significantly different by using a binomial test, which has the following characteristics:
\[H_0: \pi = \pi_0\] In our case the theoretical π will equal to 0.0909091, while the probability found in the observations is 0.2966667. We want to test if waht it observed is different from the expected value.
(bintest.poiss <- binom.test(poiss.violations, 300, p = 1-p))
##
## Exact binomial test
##
## data: poiss.violations and 300
## number of successes = 89, number of trials = 300, p-value < 2.2e-16
## alternative hypothesis: true probability of success is not equal to 0.09090909
## 95 percent confidence interval:
## 0.2455466 0.3518550
## sample estimates:
## probability of success
## 0.2966667
We can see that for the poisson, the p-value is really low 4.138435310^{-24}. This means that we will reject the null hypothesis and hence that the distribution of the violations is not following a bernoulli distributions. There has been an underestimation of the quantile we used to calculate the VaR and the values at risk, as we had quite a lot violations compared to the theoretical one. We believe that by increasing the probability in the poisson distribution we would get to a more optimal result.
We will split the data in two subgroups, the first one being the training set which we will use to estimate the VaR will comprehend the observations up to the last 300 days of the data, while the last 300 observations will be the test set to which we will compare the VaR we will find.
salesTest <- sales$Sales[(nrow(sales)-299):(nrow(sales))]
Then we will calculate the VaR for the test set using a moving window of 365 days.
In order to do so, we start by creating a function to calculate the VaR for a given window, which works in the following way: - here we use the 73% quantile to determine the threshold as a higher one would mean too few observations above the line to be able to fit a gpd and it would create too few violations, while a lower one would create the opposite problems, - window of 365 days before the observations that we want to estimate, - calculate the weights for the observations in this window which are above the threshold by taking the differences between the sales and the thershold, - fit a GPD on the weights and estimate the parameters for the shape and the scale, - same alpha as before, which is the p value we calculated, - estimation of Fbar ny dividing the number of observations above the threshold in the window and the total number of observation in the window, - put all the estimators in the formula to calculate the value at risk.
VaRest <- function(df){
#threshold
u <- 0.73 * max(df$Sales)
#Weights for the observations above the threshold chosen
W <- df %>% dplyr::filter(Sales > u) %>% dplyr::mutate(W = Sales - u) %>% dplyr::select(W)
#fitting a GPD, method amle as we want the shape parameter positive
gpd <- gPdtest::gpd.fit(W$W, method = "amle")
#extracting the parameters
#shape
beta <- gpd[1,]
#scale
xi <- gpd[2,]
#setting the alpha for the VaR to the level we have calculated
alpha <- p
#calculating Fbar with the approximation
Fbar <- nrow(W) / 300
#calculating the values at risk
(u + (beta/xi) * ( ( ( (1 - alpha) / (Fbar) ) ^ (-xi) ) - 1) )
}
Then we move on with the calculation of the value at risk for each observation in the test set. In order to do it, we will first create a vector to store the values, which will have a length of 300 as the number of observations in the test set. Then, we create a loop for the 300 observations, we want to create the window considered to calculate the value at risk by including the 365 observations that precedes the one for which we want to estimate the value at risk. Eventually, we get the different values at risk for each one of the observations.
#create a vector to store the values
VAR <- as.numeric(1:300)
for(i in 1:300){
#select window
wind <- sales %>% dplyr::slice((1913-i):(1913-i-365))
#calculate VaR
VAR[i] <- VaRest(wind)
}
VAR
## [1] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
## [9] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
## [17] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
## [25] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
## [33] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
## [41] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
## [49] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
## [57] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504
## [65] 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6504 138.6438
## [73] 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438
## [81] 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438
## [89] 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438
## [97] 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438
## [105] 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438
## [113] 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438 138.6438
## [121] 138.6438 138.6504 135.7895 135.7767 135.7895 135.7895 135.7767 135.7767
## [129] 135.7767 135.7767 135.7767 135.7885 135.8013 135.8013 135.8013 135.8013
## [137] 135.8013 135.8013 135.8013 135.8013 135.8013 135.7886 135.7886 135.7886
## [145] 135.7886 135.7886 135.7886 135.7886 135.8015 135.8015 135.8015 135.8015
## [153] 135.8015 135.8015 135.8015 135.8015 135.8015 135.8015 135.8015 135.8015
## [161] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159
## [169] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159
## [177] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159
## [185] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159
## [193] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159
## [201] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159
## [209] 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8159 135.8310
## [217] 135.8310 135.8310 135.8310 135.8310 135.8310 135.8310 135.8310 135.8481
## [225] 135.8481 135.8481 135.8481 135.8481 135.8481 135.8481 135.8667 135.8667
## [233] 135.8667 135.8667 135.8667 135.8893 135.8893 135.8893 135.8893 135.8893
## [241] 135.8893 135.8893 135.8893 135.9120 135.9120 133.0342 121.1951 121.1566
## [249] 121.1566 121.1566 121.1566 121.1951 121.2387 121.5214 121.4619 121.5214
## [257] 121.5830 121.6408 150.9121 152.5583 152.4983 152.4983 152.4983 152.4983
## [265] 152.4953 152.4927 152.4927 152.4927 152.4927 152.4927 152.4927 152.4927
## [273] 152.4927 152.4904 152.4904 152.4904 152.4904 152.4904 152.4904 152.4904
## [281] 152.4904 152.4927 152.4927 152.4927 152.4927 152.4927 152.4904 152.4904
## [289] 152.4904 152.4904 152.4904 152.4904 152.4904 152.4904 152.4904 152.4927
## [297] 152.4927 152.4927 152.4927 152.4927
We can see that there are quite a lot of repetitions in the estimation of the values at risk, this is given by the fact that the the observations contained in the windows of adjacent dates are all the same but one, as we discard the oldest one to replace it with the real value of the observation of the day before the one we are considering. Hence, if none of these two values is describing a violation (meaning that it is higher than the threshold that we have chosen), the computation of the value at risk won’t be affected, as all the paramters that we use to estimate the VaR will remain unchanged: the weights will be the same (since we consider only the values that are above the threshold) and so will be the GPD and its parameters’ estimations, the fbar is also the same as Nu stays the same and the alpha is fixed.
Now, we compare the VaR we estimated to the real values, to see how many violations there are. The should be following a bernoulli distribution with n = 300 and p = 1 - alpha.
#number of observed violations
violations <- sum(ifelse(VAR > salesTest, 0, 1))
#number of theoretical violations, given by the mean of the bernoully distribution
ev <- 300*(1 - p)
We find out that there are 21 violations.
We can compare this value to the theoretical one that we expect to get, which is 27.2727273.
We can test whether they are statistically significantly different by using a binomial test, which has the following characteristics:
\[H_0: \pi = \pi_0\] In our case the theoretical π will equal to 0.0909091, while the probability found in the observations is 0.07. We want to test if waht it observed is different from the expected value.
(bintest <- binom.test(violations, 300, p = 1-p))
##
## Exact binomial test
##
## data: violations and 300
## number of successes = 21, number of trials = 300, p-value = 0.2287
## alternative hypothesis: true probability of success is not equal to 0.09090909
## 95 percent confidence interval:
## 0.04385034 0.10501400
## sample estimates:
## probability of success
## 0.07
We can see that the p-value is quite high (0.2286663), more specifically is higher than the significance level alpha = 5%, hence we cannot reject the null hypothesis of the two proportions being the same and hence that the distribution of the violations in the windows considered follows a binomial(300, 1-p).
In this third practical, we expand the retail sales dataset of Practical 2 to include multiple products. Data for Product 1 (from Practical 2) and an additional 9 products are given in two forms: \(sales10.csv\) and \(sales10adjusted.csv\), with the same format. As the name indicates, sales volumes in the latter file have been adjusted for trend and seasonality (removing growth and cyclical patterns).
#data
sales10 <- read.csv(here::here("sales_10.csv"))
sales10adj <- read.csv(here::here("sales_10_adjusted.csv"))
Explain the challenges involved when analysing the risk of stock-out across several products. What is the danger of considering each series independently?
The aplication of an independently analysis (EVT) could cause a loss of information and therefore a lack of preparation in face of a risky event. For example, the analysis of stock out from 3 products independently are not consider risky, nevertheless, the stock out of the 3 together could be consider as a very high risky event due to the high probability of loosing a customer.
MEVT has as a main objective to find the probabilities that extreme events happen together.
The main problem would be given by the fact that the different products may not be independent even if we consider them to be, it can be seen if one of them has an extreme value in the demand and so does another one, every time at the same date. This can happen as the products may be related, as for example could be bread and butter, meaning that if there is an increase in the demand of one product, so there will be in the demand for the other one. Another possibility is the fact that, especially for extreme values, so an abrupt change in the demand, could be determined by external factors, such as a natural disaster, that would make the demand for all sort of products increase from one day to the other. Indeed, it could be a problem, as by analysing the stock-outs, which are extreme values, independently one from the other, we could miss some important relationship that are between them and not being able to assess the risk in the best way possible.
Propose one graphical and one numerical method of detecting dependence of extreme values of the demand across several products. Apply your chosen methods to the adjusted sales data and identify groups of related products (if any exist).
| n° | Analysis | Method |
|---|---|---|
| 1 | Graphical | We will apply two methods: 1) A plot with all the products distributions over time, in which we will set a threshold as an arbitrary value, in order to enable us the study the ones which are above it and see if there is repeated pattern (for example having always two products that are over the threshold in two close dates) 2) A chi-plot among all the different couples of products, so that we can study the distribution of the tails of the products’ demands. |
| 2 | Numerical | We will compute the chibar values, which will confirm us the tail dependences among all the different couples of products. |
| 3 | Aggregation | In this step we will focus in the grouping of products, in order to do so, we will use the following 2 method: 1) We will calculate the mean and variance for each product and display them on a graph, to see if they have any other product close to them 2) Then, we will compute the Kendall-Taus distances among the products and see the clustering of the products based on it. |
First, let’s create a plot with all the products’ demands distibutions.
## integer(0)
In order to determine the threshold for the extreme values, we will take the 98% quantile of the distribution of the sales of each product, as we could find very different distributions and taking only one threshold overall would not give us any important information.
As we can see, we have some observations in each product that are above the line. We are particularly interested in seeing if there is a repetition in the observations that are above the lines that happen around the same time for different products.
Now, we will determine the number of the observation in which we pass the line for each products and we try to see if there are some days that are close enough or repetead.
| Product 1 | Product 2 | Product 3 | Product 4 | Product 5 | Product 6 | Product 7 | Product 8 | Product 9 | Product 10 |
|---|---|---|---|---|---|---|---|---|---|
| 240 | 300 | 205 | 66 | 299 | 8 | 209 | 242 | 102 | 39 |
| 299 | 301 | 301 | 154 | 330 | 9 | 223 | 253 | 104 | 53 |
| 329 | 329 | 329 | 301 | 630 | 325 | 325 | 321 | 226 | 101 |
| 336 | 332 | 332 | 332 | 663 | 330 | 349 | 334 | 249 | 193 |
| 337 | 382 | 382 | 333 | 696 | 337 | 365 | 345 | 282 | 226 |
| 664 | 387 | 387 | 336 | 940 | 522 | 380 | 372 | 283 | 248 |
| 696 | 665 | 388 | 584 | 946 | 664 | 413 | 413 | 307 | 255 |
| 703 | 696 | 576 | 665 | 947 | 696 | 468 | 436 | 333 | 273 |
| 938 | 888 | 664 | 694 | 976 | 703 | 471 | 611 | 437 | 282 |
| 960 | 933 | 696 | 696 | 977 | 945 | 497 | 618 | 499 | 283 |
| 962 | 1061 | 1061 | 699 | 1039 | 955 | 584 | 788 | 522 | 285 |
| 1018 | 1226 | 1281 | 888 | 1042 | 1038 | 701 | 974 | 603 | 317 |
| 1019 | 1274 | 1304 | 1060 | 1044 | 1042 | 758 | 1263 | 736 | 432 |
| 1039 | 1278 | 1426 | 1112 | 1063 | 1043 | 771 | 1264 | 744 | 470 |
| 1041 | 1279 | 1428 | 1252 | 1068 | 1044 | 779 | 1335 | 745 | 590 |
| 1044 | 1312 | 1526 | 1253 | 1085 | 1068 | 784 | 1339 | 829 | 614 |
| 1068 | 1399 | 1554 | 1275 | 1398 | 1136 | 828 | 1690 | 1166 | 645 |
| 1136 | 1426 | 1790 | 1399 | 1402 | 1191 | 918 | 1728 | 1261 | 683 |
| 1664 | 1428 | 1791 | 1425 | 1475 | 1526 | 1435 | 1736 | 1382 | 1166 |
| 1762 | 1791 | 1793 | 1428 | 1762 | 1763 | 1517 | 1910 | 1746 | 1382 |
We can see that especially for the dates of products 2 and 3 we have some repetitions of the same values, also for products 8 and 7 it looks like there is some similitudes or closeness of the values.
The other method we can use is to create the chi plot for all the couples of products to assess the distributions of their tails and see if there is any correlation.
Key takeaway
We can see that in the highest quantiles we have a straight line for the majority of products, only for a couple of them (like p1p3 and p5p6), there seems to be a line that is not centered around zero, but to a higher value, meaning that there could be some correlation among them, which could be a problem.
While for the mathematical method, we will construct a matrix containing the chi bar values for all the couples of products.
| p2 | p3 | p4 | p5 | p6 | p7 | p8 | p9 | p10 | |
|---|---|---|---|---|---|---|---|---|---|
| p1 | 0.218 | 0.181 | 0.200 | 0.442 | 0.442 | 0.039 | 0.218 | 0.094 | -0.072 |
| p2 | 0.000 | 0.508 | 0.430 | 0.236 | 0.140 | 0.161 | 0.068 | 0.094 | 0.007 |
| p3 | 0.000 | 0.000 | 0.406 | 0.218 | 0.218 | 0.118 | 0.094 | 0.007 | 0.068 |
| p4 | 0.000 | 0.000 | 0.000 | 0.161 | 0.140 | 0.039 | 0.161 | 0.094 | -0.072 |
| p5 | 0.000 | 0.000 | 0.000 | 0.000 | 0.430 | 0.039 | 0.118 | 0.039 | -0.029 |
| p6 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | -0.029 | -0.029 | 0.007 | -0.127 |
| p7 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.094 | 0.039 | 0.007 |
| p8 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.068 | 0.118 |
| p9 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.328 |
We can see that the highest values appear to be for products 2 with product 3 and 4, product 1 with product 5 and 6, product 3 with 4, products 5 with 6. This is very much in line with we have found in the chi plot that we created above.
To investigate a possible clustering of the different products, we will create create a graph with the mean and variance for each product.
Key takeaway What is interesting to notice, is that the clustering are the same if we consider the whole distribution and only the most extreme values. Netherless, the correlation are stronger in the extreme values, thus we confirm that an independent analysis could potentialy lose important information.
In both cases, we can see that there is some similarity in some products, and some differences in others, especially product 1 appreas to be very different from the others, while products 9, 10 7, 8, 4, 3 and 2 appear very similar (the cluster in green). However, the number of clusters that should be used is not very clear, as the decrease of total within sum of sqaure is quite linear.
In view of your answer to (b), would you apply a Gaussian copula model to these data?
As well as the EVT the Gaussian model does not caputre all the information for the extreme values. Keeping this concept for multivariate extreme value theory we can see that in this case the corelation is higher, thus we cannot apply this model to fit the extreme value of the products’ distribution.
As we could find some depence in the distribution of the tails of our products’ demand, we believe that a gaussian distribution would not be optimal, as it is not able to consider the depences in the tails of the distribution.
Fit a multivariate model to the adjusted sales data and estimate the 95% Demand at Risk for the sum:
\[S = X^{(1)} + X^{(2)} + ... + X^{(10)}\]
- (Hard) Fit a multivariate model to the adjusted sales data and estimate the 95% Demand at Risk for the sum: \[S = X^{(1)} + X^{(2)} + ... + X^{(10)}\] You can do so with the following steps:
# These functions can be used to transform data to uniform scale,
# when the data are analysed separately below and above a threshold
# e.g. applying the GPD to values above a high threshold and treating
# those below the threshold with the empirical CDF.
transform_to_uniform <- function(y,
cdf,
threshold,
...) {
# ...: arguments to `cdf`
above <- y > threshold
# proportion of exceedances
p <- mean(above)
# under the threshold, use the empirical CDF
# see ?ecdf
ecdf_below <- ecdf(y[!above])
empirical <- ecdf_below(y[!above])
# above, apply cdf
theoretical <- cdf(y[above], mean.sales10adj)
transformed <- numeric(length = length(y))
# "Glue together" the empirical and theoretical parts by rescaling them:
# empirical e.g. N(0,1)
# [-----------------------|-------------]
# 0 1-p 1
transformed[!above] <- (1 - p) * empirical # (1)
transformed[above] <- 1 - p + p * theoretical # (2)
return(list(
transformed = transformed,
ecdf = ecdf_below,
prop = p
))
}
inverse_transform <- function(u, ecdf, quantile_function, p, ...) {
# ...: arguments to quantile_function
# Basically, undo the transformations (1-2) above
above <- u > 1 - p
original_scale <- numeric(length = length(u))
original_scale[!above] <- quantile(ecdf, u[!above] / (1 - p))
original_scale[above] <-
quantile_function((u[above] - (1 - p)) / p, ...)
return(original_scale)
}
Transform the sales for each product to a uniform scale, that is, compute \[F_i(X^{(i)})\] for each X in 1, …, n being each of the product sales and F being the estimated Extreme Value distributions.
Fit a multivariate copula of your choice using the copula::fitCopula function.
library(copula)
## Warning: package 'copula' was built under R version 3.6.2
unif <- cbind(unif1$transformed, unif2$transformed, unif3$transformed, unif4$transformed, unif5$transformed, unif6$transformed, unif7$transformed, unif8$transformed, unif9$transformed, unif10$transformed)
gum.fit <- gumbelCopula(5, dim=10)
clay.fit <- claytonCopula(4, dim=10)
fit.gum <- fitCopula(gum.fit, method = "itau", data=unif)
## Warning in .local(copula, tau, ...): For the Gumbel copula, tau must be >= 0.
## Replacing negative values by 0.
summary(fit.gum)
## Call: fitCopula(copula, data = data, method = "itau")
## Fit based on "inversion of Kendall's tau" and 1913 10-dimensional observations.
## Gumbel copula, dim. d = 10
## Estimate Std. Error
## alpha 1.04 NA
Simulate from your fitted copula and transform the values back to their original scales (i.e. undo the transformation in i.)
Compute the simulated values of S and their 95% Value at Risk.
## 'data.frame': 50 obs. of 1 variable:
## $ rev10: num 1.85 3.67 2.38 1.38 1.5 ...
| rev1 | rev2 | rev3 | rev4 | rev5 | rev6 | rev7 | rev8 | rev9 | rev10 |
|---|---|---|---|---|---|---|---|---|---|
| 1.833397 | 1.721064 | 1.674045 | 2.215008 | 1.289980 | 1.6627326 | 1.729129 | 1.754176 | 1.383327 | 1.848719 |
| 1.389896 | 2.381227 | 2.229203 | 1.342169 | 1.162100 | 1.1612029 | 2.654195 | 2.021832 | 1.939627 | 3.669373 |
| 1.613484 | 1.321973 | 1.447119 | 1.644108 | 1.619078 | 1.1096222 | 1.435127 | 1.564317 | 1.424517 | 2.384326 |
| 1.785663 | 1.704681 | 1.697891 | 1.674859 | 1.839798 | 2.1217804 | 1.311406 | 1.458367 | 1.444699 | 1.379314 |
| 1.208821 | 1.198224 | 1.581209 | 1.519660 | 1.101783 | 1.4266296 | 1.190012 | 1.490266 | 1.510652 | 1.496257 |
| 1.634306 | 1.590843 | 1.607376 | 1.675039 | 1.295762 | 1.7613302 | 3.214484 | 1.391533 | 1.703005 | 1.556532 |
| 1.533151 | 1.605448 | 1.601066 | 1.915172 | 1.079881 | 1.2236936 | 1.330645 | 1.609578 | 2.456184 | 1.364795 |
| 1.190154 | 2.004760 | 1.755000 | 1.918430 | 1.271633 | 1.2601710 | 1.193537 | 1.671078 | 1.876376 | 1.357561 |
| 1.233841 | 1.562571 | 1.978563 | 2.452361 | 1.281525 | 1.0810867 | 2.127564 | 1.349537 | 1.920688 | 1.536425 |
| 1.496956 | 1.982314 | 2.487770 | 1.444414 | 2.045698 | 1.2799271 | 2.037815 | 1.575224 | 1.249127 | 1.366893 |
| 2.113782 | 1.321333 | 1.474030 | 1.897923 | 1.423585 | 1.2549105 | 2.865743 | 1.721762 | 1.800719 | 1.202008 |
| 1.607884 | 2.760839 | 1.988598 | 1.933120 | 1.177907 | 2.4900720 | 1.756004 | 1.597384 | 1.611118 | 1.663281 |
| 2.873292 | 1.842880 | 1.519386 | 2.084405 | 1.643726 | 1.2696415 | 1.459083 | 1.387173 | 1.477006 | 1.200218 |
| 1.210557 | 1.101677 | 1.461695 | 1.356342 | 1.226587 | 2.0802226 | 1.816353 | 1.888337 | 1.265161 | 1.740980 |
| 1.585486 | 3.635026 | 1.416224 | 1.432849 | 1.169583 | 2.0454140 | 1.214606 | 1.546852 | 1.169102 | 1.778391 |
| 1.331183 | 1.123347 | 1.492611 | 2.248024 | 1.785575 | 1.7228191 | 1.311531 | 1.902817 | 1.256973 | 1.841823 |
| 1.256543 | 1.688945 | 1.376746 | 1.544800 | 1.525465 | 0.9615920 | 1.199900 | 1.537253 | 1.552004 | 1.660870 |
| 1.296168 | 1.527678 | 1.485815 | 1.557216 | 1.664586 | 1.9478886 | 2.309075 | 1.417349 | 1.952647 | 1.339120 |
| 1.181620 | 2.049961 | 2.321236 | 1.355867 | 1.337488 | 2.2613678 | 1.656007 | 1.377070 | 1.485760 | 1.945392 |
| 2.100210 | 1.077789 | 1.433739 | 1.592683 | 1.418845 | 1.9574048 | 1.806083 | 1.486208 | 2.004329 | 1.487589 |
| 2.061404 | 1.913821 | 1.526011 | 1.564677 | 1.936295 | 1.6178485 | 1.749009 | 2.314871 | 2.745721 | 1.265690 |
| 1.517808 | 2.310673 | 1.970856 | 1.690513 | 2.531249 | 1.1877103 | 1.629864 | 1.716186 | 1.393269 | 1.860563 |
| 1.733879 | 2.434929 | 2.002785 | 1.615377 | 1.387149 | 1.9369006 | 2.445183 | 1.914617 | 1.360017 | 2.552996 |
| 1.089229 | 1.686079 | 1.450542 | 2.046041 | 2.275769 | 1.7553744 | 1.859622 | 1.632278 | 1.644185 | 1.424989 |
| 1.442043 | 1.457383 | 1.398268 | 1.836652 | 1.536172 | 0.9777846 | 1.261231 | 1.851284 | 1.188082 | 2.038215 |
| 1.659330 | 1.776774 | 1.362589 | 1.640176 | 1.108892 | 1.2932010 | 1.259842 | 2.013255 | 1.482082 | 1.527087 |
| 2.070678 | 1.731888 | 2.068902 | 1.868487 | 1.719934 | 1.8460408 | 1.547536 | 2.343498 | 2.364856 | 1.509815 |
| 1.474110 | 1.288580 | 1.660731 | 1.388430 | 2.238346 | 1.0857424 | 1.239393 | 1.507835 | 1.746296 | 1.285536 |
| 1.352808 | 1.179803 | 2.185961 | 1.485179 | 1.554753 | 1.1317846 | 2.322521 | 1.376496 | 2.409522 | 1.659858 |
| 2.278935 | 1.112572 | 1.419617 | 1.421075 | 1.114297 | 1.6818651 | 1.509684 | 1.583086 | 1.365620 | 1.382108 |
| 1.370466 | 2.033608 | 1.692362 | 1.491893 | 1.742530 | 1.4601652 | 1.220284 | 1.840711 | 1.532514 | 1.341479 |
| 2.271847 | 1.241626 | 2.203395 | 1.500229 | 1.344072 | 0.9984627 | 1.136363 | 1.975560 | 1.423543 | 1.359376 |
| 1.275806 | 1.911564 | 2.329551 | 1.459843 | 1.085427 | 1.3701531 | 1.320556 | 2.467837 | 2.056330 | 1.359570 |
| 2.040185 | 1.360157 | 1.344735 | 2.155631 | 1.550392 | 1.1173153 | 1.438916 | 1.965333 | 1.369547 | 2.068387 |
| 1.358254 | 1.464846 | 2.032040 | 1.799191 | 1.976331 | 1.2252450 | 1.448942 | 1.272767 | 1.345834 | 2.668137 |
| 1.214027 | 1.105772 | 1.565685 | 2.043336 | 1.645061 | 1.8596732 | 1.326358 | 1.691640 | 1.362605 | 2.307718 |
| 1.831860 | 2.649468 | 1.902648 | 1.695757 | 2.128824 | 2.8308023 | 1.248845 | 1.895207 | 1.390897 | 1.223339 |
| 1.946117 | 1.639319 | 1.518642 | 1.798334 | 1.470069 | 2.4631674 | 1.379029 | 2.494113 | 1.425435 | 1.330472 |
| 1.647007 | 1.389182 | 1.393164 | 1.372668 | 1.168176 | 1.4802188 | 1.445032 | 1.802889 | 1.339266 | 1.513919 |
| 2.111566 | 1.621631 | 2.088837 | 1.377249 | 2.392124 | 1.0427114 | 3.191552 | 1.875108 | 2.859082 | 1.236406 |
| 2.000295 | 1.516690 | 1.353543 | 2.443163 | 1.752450 | 3.0232093 | 1.307114 | 1.446175 | 1.901184 | 1.478862 |
| 1.430779 | 1.093592 | 1.640822 | 1.374975 | 1.842401 | 1.7282744 | 1.288309 | 1.915499 | 1.637622 | 1.528434 |
| 2.191231 | 1.139949 | 2.386237 | 2.314526 | 1.615161 | 1.1275522 | 1.909695 | 2.169519 | 1.483068 | 1.271424 |
| 2.358958 | 1.933521 | 2.134471 | 1.374123 | 1.654963 | 2.2836139 | 1.991499 | 1.647774 | 2.126623 | 1.221293 |
| 1.525797 | 2.363068 | 1.526451 | 1.380267 | 2.215259 | 1.9524406 | 1.788769 | 1.877078 | 1.670259 | 1.261554 |
| 1.901465 | 1.395117 | 1.653835 | 1.401864 | 1.668551 | 2.5688082 | 1.586492 | 1.833891 | 2.193103 | 1.581098 |
| 1.454989 | 1.165118 | 1.545595 | 2.484355 | 2.508851 | 1.2133469 | 1.145841 | 1.721651 | 1.805971 | 1.391955 |
| 1.721983 | 1.141894 | 3.133116 | 1.607453 | 1.309691 | 0.9835133 | 1.310715 | 2.029665 | 1.703807 | 2.299050 |
| 2.038518 | 1.195019 | 2.221920 | 1.631545 | 1.242957 | 1.3963297 | 1.266906 | 1.292525 | 1.666654 | 1.737291 |
| 1.322076 | 1.525024 | 1.427202 | 1.531989 | 1.465071 | 1.4171284 | 1.581289 | 1.596690 | 1.173473 | 1.476023 |
| rev1 | rev2 | rev3 | rev4 | rev5 | rev6 | rev7 | rev8 | rev9 | rev10 |
|---|---|---|---|---|---|---|---|---|---|
| 1.833397 | 1.721064 | 1.674045 | 2.215008 | 1.289980 | 1.6627326 | 1.729129 | 1.754176 | 1.383327 | 1.848719 |
| 1.389896 | 2.381227 | 2.229203 | 1.342169 | 1.162100 | 1.1612029 | 2.654195 | 2.021832 | 1.939627 | 3.669373 |
| 1.613484 | 1.321973 | 1.447119 | 1.644108 | 1.619078 | 1.1096222 | 1.435127 | 1.564317 | 1.424517 | 2.384326 |
| 1.785663 | 1.704681 | 1.697891 | 1.674859 | 1.839798 | 2.1217804 | 1.311406 | 1.458367 | 1.444699 | 1.379314 |
| 1.208821 | 1.198224 | 1.581209 | 1.519660 | 1.101783 | 1.4266296 | 1.190012 | 1.490266 | 1.510652 | 1.496257 |
| 1.634306 | 1.590843 | 1.607376 | 1.675039 | 1.295762 | 1.7613302 | 3.214484 | 1.391533 | 1.703005 | 1.556532 |
| 1.533151 | 1.605448 | 1.601066 | 1.915172 | 1.079881 | 1.2236936 | 1.330645 | 1.609578 | 2.456184 | 1.364795 |
| 1.190154 | 2.004760 | 1.755000 | 1.918430 | 1.271633 | 1.2601710 | 1.193537 | 1.671078 | 1.876376 | 1.357561 |
| 1.233841 | 1.562571 | 1.978563 | 2.452361 | 1.281525 | 1.0810867 | 2.127564 | 1.349537 | 1.920688 | 1.536425 |
| 1.496956 | 1.982314 | 2.487770 | 1.444414 | 2.045698 | 1.2799271 | 2.037815 | 1.575224 | 1.249127 | 1.366893 |
| 2.113782 | 1.321333 | 1.474030 | 1.897923 | 1.423585 | 1.2549105 | 2.865743 | 1.721762 | 1.800719 | 1.202008 |
| 1.607884 | 2.760839 | 1.988598 | 1.933120 | 1.177907 | 2.4900720 | 1.756004 | 1.597384 | 1.611118 | 1.663281 |
| 2.873292 | 1.842880 | 1.519386 | 2.084405 | 1.643726 | 1.2696415 | 1.459083 | 1.387173 | 1.477006 | 1.200218 |
| 1.210557 | 1.101677 | 1.461695 | 1.356342 | 1.226587 | 2.0802226 | 1.816353 | 1.888337 | 1.265161 | 1.740980 |
| 1.585486 | 3.635026 | 1.416224 | 1.432849 | 1.169583 | 2.0454140 | 1.214606 | 1.546852 | 1.169102 | 1.778391 |
| 1.331183 | 1.123347 | 1.492611 | 2.248024 | 1.785575 | 1.7228191 | 1.311531 | 1.902817 | 1.256973 | 1.841823 |
| 1.256543 | 1.688945 | 1.376746 | 1.544800 | 1.525465 | 0.9615920 | 1.199900 | 1.537253 | 1.552004 | 1.660870 |
| 1.296168 | 1.527678 | 1.485815 | 1.557216 | 1.664586 | 1.9478886 | 2.309075 | 1.417349 | 1.952647 | 1.339120 |
| 1.181620 | 2.049961 | 2.321236 | 1.355867 | 1.337488 | 2.2613678 | 1.656007 | 1.377070 | 1.485760 | 1.945392 |
| 2.100210 | 1.077789 | 1.433739 | 1.592683 | 1.418845 | 1.9574048 | 1.806083 | 1.486208 | 2.004329 | 1.487589 |
| 2.061404 | 1.913821 | 1.526011 | 1.564677 | 1.936295 | 1.6178485 | 1.749009 | 2.314871 | 2.745721 | 1.265690 |
| 1.517808 | 2.310673 | 1.970856 | 1.690513 | 2.531249 | 1.1877103 | 1.629864 | 1.716186 | 1.393269 | 1.860563 |
| 1.733879 | 2.434929 | 2.002785 | 1.615377 | 1.387149 | 1.9369006 | 2.445183 | 1.914617 | 1.360017 | 2.552996 |
| 1.089229 | 1.686079 | 1.450542 | 2.046041 | 2.275769 | 1.7553744 | 1.859622 | 1.632278 | 1.644185 | 1.424989 |
| 1.442043 | 1.457383 | 1.398268 | 1.836652 | 1.536172 | 0.9777846 | 1.261231 | 1.851284 | 1.188082 | 2.038215 |
| 1.659330 | 1.776774 | 1.362589 | 1.640176 | 1.108892 | 1.2932010 | 1.259842 | 2.013255 | 1.482082 | 1.527087 |
| 2.070678 | 1.731888 | 2.068902 | 1.868487 | 1.719934 | 1.8460408 | 1.547536 | 2.343498 | 2.364856 | 1.509815 |
| 1.474110 | 1.288580 | 1.660731 | 1.388430 | 2.238346 | 1.0857424 | 1.239393 | 1.507835 | 1.746296 | 1.285536 |
| 1.352808 | 1.179803 | 2.185961 | 1.485179 | 1.554753 | 1.1317846 | 2.322521 | 1.376496 | 2.409522 | 1.659858 |
| 2.278935 | 1.112572 | 1.419617 | 1.421075 | 1.114297 | 1.6818651 | 1.509684 | 1.583086 | 1.365620 | 1.382108 |
| 1.370466 | 2.033608 | 1.692362 | 1.491893 | 1.742530 | 1.4601652 | 1.220284 | 1.840711 | 1.532514 | 1.341479 |
| 2.271847 | 1.241626 | 2.203395 | 1.500229 | 1.344072 | 0.9984627 | 1.136363 | 1.975560 | 1.423543 | 1.359376 |
| 1.275806 | 1.911564 | 2.329551 | 1.459843 | 1.085427 | 1.3701531 | 1.320556 | 2.467837 | 2.056330 | 1.359570 |
| 2.040185 | 1.360157 | 1.344735 | 2.155631 | 1.550392 | 1.1173153 | 1.438916 | 1.965333 | 1.369547 | 2.068387 |
| 1.358254 | 1.464846 | 2.032040 | 1.799191 | 1.976331 | 1.2252450 | 1.448942 | 1.272767 | 1.345834 | 2.668137 |
| 1.214027 | 1.105772 | 1.565685 | 2.043336 | 1.645061 | 1.8596732 | 1.326358 | 1.691640 | 1.362605 | 2.307718 |
| 1.831860 | 2.649468 | 1.902648 | 1.695757 | 2.128824 | 2.8308023 | 1.248845 | 1.895207 | 1.390897 | 1.223339 |
| 1.946117 | 1.639319 | 1.518642 | 1.798334 | 1.470069 | 2.4631674 | 1.379029 | 2.494113 | 1.425435 | 1.330472 |
| 1.647007 | 1.389182 | 1.393164 | 1.372668 | 1.168176 | 1.4802188 | 1.445032 | 1.802889 | 1.339266 | 1.513919 |
| 2.111566 | 1.621631 | 2.088837 | 1.377249 | 2.392124 | 1.0427114 | 3.191552 | 1.875108 | 2.859082 | 1.236406 |
| 2.000295 | 1.516690 | 1.353543 | 2.443163 | 1.752450 | 3.0232093 | 1.307114 | 1.446175 | 1.901184 | 1.478862 |
| 1.430779 | 1.093592 | 1.640822 | 1.374975 | 1.842401 | 1.7282744 | 1.288309 | 1.915499 | 1.637622 | 1.528434 |
| 2.191231 | 1.139949 | 2.386237 | 2.314526 | 1.615161 | 1.1275522 | 1.909695 | 2.169519 | 1.483068 | 1.271424 |
| 2.358958 | 1.933521 | 2.134471 | 1.374123 | 1.654963 | 2.2836139 | 1.991499 | 1.647774 | 2.126623 | 1.221293 |
| 1.525797 | 2.363068 | 1.526451 | 1.380267 | 2.215259 | 1.9524406 | 1.788769 | 1.877078 | 1.670259 | 1.261554 |
| 1.901465 | 1.395117 | 1.653835 | 1.401864 | 1.668551 | 2.5688082 | 1.586492 | 1.833891 | 2.193103 | 1.581098 |
| 1.454989 | 1.165118 | 1.545595 | 2.484355 | 2.508851 | 1.2133469 | 1.145841 | 1.721651 | 1.805971 | 1.391955 |
| 1.721983 | 1.141894 | 3.133116 | 1.607453 | 1.309691 | 0.9835133 | 1.310715 | 2.029665 | 1.703807 | 2.299050 |
| 2.038518 | 1.195019 | 2.221920 | 1.631545 | 1.242957 | 1.3963297 | 1.266906 | 1.292525 | 1.666654 | 1.737291 |
| 1.322076 | 1.525024 | 1.427202 | 1.531989 | 1.465071 | 1.4171284 | 1.581289 | 1.596690 | 1.173473 | 1.476023 |